/*
 * Rho.cpp
 *
 *  Created on: May 24, 2011
 *      Author: marchi
 */

#include "Rho.h"
#include <iostream>
#include <cmath>


const rvec Rho::cube[VERTEX]={{0,0,0},{1,0,0},{0,1,0},{0,0,1},{1,1,0},{1,0,1},{0,1,1},{1,1,1}};
Rho & Rho::operator=(const Grid & x){
	Grid & My=*this;
	My=x;
	return *this;
}
Rho & Rho::operator=(const double & x){
	Grid & My=*this;
	My=x;
	return *this;
}
void Rho::TriLinear(real gx, real gy, real gz, real q){
	int n;
	real dd;
	real sum=0;
	for(n=0;n<VERTEX;n++){
		dd=TriDiff(gx,cube[n][XX])*TriDiff(gy,cube[n][YY])*TriDiff(gz,cube[n][ZZ]);
		qq[n]=dd*q;
		sum+=qq[n];
	}
}
void Rho::Density(Charges & x){
// int get_charges(real*** Rho_slt, real*** Rho_sol,real* Charges,int nx,int ny,int nz,int nind, atom_id *cindex, t_trxframe* fr){

	try{
		if(!x.test()) throw "Cannot compute density if charges not initialized.";
	}
	catch(char *s){
		std::cout << s << std::endl;
	}
	int nx0=static_cast<int>(nnx);
	int ny0=static_cast<int>(nny);
	int nz0=static_cast<int>(nnz);

	int natoms=x.getNR();
	array2<real> xa(natoms,DIM);

	real x1,y1,z1,r1,s1,t1,gx,gy,gz;
	int mx,my,mz;
	int ox,oy,oz,tx,ty,tz;
	double DV=Rho::getDV();
	*this=0.0;
	for (int i = 0; i < natoms; i++) {
		xa[i][XX] = oc[XX][XX] * x[i][XX] + oc[XX][YY] * x[i][YY] + oc[XX][ZZ] * x[i][ZZ];
		xa[i][YY] = oc[YY][XX] * x[i][XX] + oc[YY][YY] * x[i][YY] + oc[YY][ZZ] * x[i][ZZ];
		xa[i][ZZ] = oc[ZZ][XX] * x[i][XX] + oc[ZZ][YY] * x[i][YY] + oc[ZZ][ZZ] * x[i][ZZ];
	}
	for (int i = 0; i < natoms; i++) {
		x1 = xa[i][XX];
		y1 = xa[i][YY];
		z1 = xa[i][ZZ];
		r1 = static_cast<real> (nx0 * (x1 - rint(x1 - 0.5)));
		s1 = static_cast<real> (ny0 * (y1 - rint(y1 - 0.5)));
		t1 = static_cast<real> (nz0 * (z1 - rint(z1 - 0.5)));
		mx = static_cast<int>(r1);
		my = static_cast<int>(s1);
		mz = static_cast<int>(t1);
		gx = r1 - (real) mx;
		gy = s1 - (real) my;
		gz = t1 - (real) mz;
		real chg=x.getQ(i);
		TriLinear(gx, gy, gz, chg);
		for (int n = 0; n < VERTEX; n++) {
			tx = static_cast<int>(cube[n][XX]);
			ty = static_cast<int>(cube[n][YY]);
			tz = static_cast<int>(cube[n][ZZ]);
			ox = (nx0 > mx + tx) ? mx + tx : mx + tx - nx0;
			oy = (ny0 > my + ty) ? my + ty : my + ty - ny0;
			oz = (nz0 > mz + tz) ? mz + tz : mz + tz - nz0;
			(*this)[ox][oy][oz] += static_cast<double>(qq[n])/DV;
		}
	}
	xa.Deallocate();
}
